---
title: "Modeling (2020 Only)"
subtitle: Interaction Figures
date: "`r format(Sys.time(), '%B %d, %Y')`"
output:
  pdf_document:
    toc: true
    toc_depth: 4
params:
  save_tex: true
---

```{r, echo=FALSE}
knitr::opts_chunk$set(echo=FALSE,
                      message=FALSE, 
                      warning=FALSE)
```

```{r}
filepath_replication_data <- '[ENTER DIRECTORY HERE]'
```

```{r lib, message=FALSE}
library(scales)
library(tidyverse)
library(lfe)
library(knitr)
library(haven)
library(broom)
library(knitr)
library(kableExtra)
library(gridExtra)
library(ggpubr)
```

```{r source_scripts}
source("helper-objects.R")
source("helper-functions.R")
```

```{r LOAD_SURVEY_DATA_dat, message=FALSE}
dat <- read_dta(file = file.path(filepath_replication_data,
                                 "final_replication_survey_dataset.dta")) %>%
  filter(wave %in% 1:4) %>% 
  rename(Q147_would_get_covid_vaccine_yes_binary   = Q147_would_get_covid_vaccine_yb,
         Q146_flu_vaccine_2020_yes_binary          = Q146_flu_vaccine_2020_yes_binar,
         zip_hh_inc_2018_median_range_middle       = zip_hh_inc_2018_median_range_mi, 
         zip_educ_2018_high_school_or_less         = zip_educ_2018_high_school_or_le, 
         demo_household_income_range_middle        = demo_household_income_range_mid,
         demo_race_ethnicity_binned_hispanic       = demo_race_ethnicity_binned_hisp,
         demo_household_income_top_quartile_binary = demo_household_income_top_quart,
         zip_hh_inc_2018_prop_top_quartile         = zip_hh_inc_2018_prop_top_quarti,
         Q10x_plan_tried_vaccinated_binary         = Q10x_plan_tried_vaccinated_bina
          ) %>% 
  mutate_if(is.labelled, as_factor)
```

```{r data_for_model_list}
data_for_model_list <- dat %>% 
  # Make long
  pivot_longer(cols = c("Q31_mask_wearing_always_binary", 
                        "Q147_would_get_covid_vaccine_yes_binary",
                        "Q146_flu_vaccine_2020_yes_binary",
                        "Q145_flu_vaccine_last_binary"),
               names_to = "dv_name", 
               values_to = "dv_value")
```

```{r dat_long_dvs}
dat_long_dvs <- dat %>% 
  select(
    weights,
    wave,
    demo_zip_scrambled,
    prop_r,
    # DVs
    Q31_mask_wearing_always_binary,
    Q147_would_get_covid_vaccine_yes_binary,
    Q146_flu_vaccine_2020_yes_binary,
    Q145_flu_vaccine_last_binary,
    # IVs
    party_with_independents
  ) %>% 
  gather(key = dv_name, value = dv_value, 
         Q31_mask_wearing_always_binary,
         Q147_would_get_covid_vaccine_yes_binary,
         Q146_flu_vaccine_2020_yes_binary,
         Q145_flu_vaccine_last_binary
         ) %>% 
    mutate(dv_name_original = dv_name,
         dv_name          = str_remove(dv_name, "^Q[0-9]*_"),
         dv_name          = str_remove(dv_name, "_binary$"),
         dv_name          = str_remove(dv_name, "^activity_")) 
```

```{r dv_name_vec_four_primary}
dv_name_vec_four_primary <- c(
  "Q31_mask_wearing_always_binary" = "Always Wear Mask",
  "Q147_would_get_covid_vaccine_yes_binary" = "COVID Vaccine",
  "Q146_flu_vaccine_2020_yes_binary" = "Flu Vaccine 2020",
  "Q145_flu_vaccine_last_binary" = "Flu Vaccine 2019"
)
```

```{r remove_objects_for_memory_reasons, eval=FALSE}
# remove objects that aren't needed. R may run out of space fast
ls_object_class <- ls() %>% 
  set_names(., .) %>% 
  map( ~ class(get(.)))

ls_object_class_df <- ls_object_class %>% 
  keep( ~ any(. == "data.frame"))

rm_keep_these_df <- c(
  "data_for_model_list"
)

rm(list = setdiff(names(ls_object_class_df), rm_keep_these_df))
res_gc <- purrr::quietly(gc)()
```


# GOP Interactions

## Run Model

```{r formula_zc_interactions_white_college_income}
formula_zc_interactions_white_college_income <- 
  dv_value ~ republican + 
             republican:prop_r + 
             republican:zip_white +
             republican:zip_educ_2018_college_and_above +
             republican:zip_hh_inc_2018_median_range_middle +
             log(deaths_14_days_county_pc_1k + 0.1) +
             # Time Fixed Effects
             as.factor(wave) +
             # Individual Controls
             factor(demo_gender) +
             demo_household_income_range_middle +
             demo_household_income_missing +
             as_factor(demo_education_binned) +
             demo_age +
             demo_race_ethnicity_binned_hispanic +
             Q7_health_dx_sum | 
             # FE + SE
             demo_zip_scrambled | 0 | demo_zip_scrambled
```

```{r model_zc_interactions_white_college_income}
model_zc_interactions_white_college_income <- data_for_model_list %>%
  # Run models
  split(.$dv_name) %>% 
  map( ~ felm(data = drop_na(., party),
              weights = drop_na(., party)$weights,
              formula = formula_zc_interactions_white_college_income))
```

## Save Tex File - All Coefficients (Table S8)

```{r model_zc_interactions_white_college_income_kable_ready}
model_zc_interactions_white_college_income_kable_ready <- 
  tidy_felm_prep_kable(
    input_model_list = model_zc_interactions_white_college_income,
    input_term_names = df_display_names_regressions,
    dv_name_vec = dv_name_vec_four_primary
  )
```

```{r model_zc_interactions_white_college_income_kable_fe_info}
model_zc_interactions_white_college_income_kable_fe_info <- create_kable_felm(
  data = model_zc_interactions_white_college_income_kable_ready,
  dv_name_vec = dv_name_vec_four_primary,
  input_term_names = df_display_names_regressions, 
  packed_rows = FALSE, 
  escape_first_space = TRUE,
  caption = "Within-ZIP Code Models of Neighborhood Partisanship (All Coefficients)",
  fixed_effect_info = list(
    vars = names(dv_name_vec_four_primary),
    fe   = c("\\ ZIP Code Fixed Effect"),
    cell_value = "\\ \\ \\ \\ \\ \\ \\checkmark"
  )
)
```

```{r SAVE_model_zc_interactions_white_college_income_kable_fe_info}
save_kable(
  x = model_zc_interactions_white_college_income_kable_fe_info,
  file = "tex-files/model_zc_interactions_white_college_income_kable_fe_info.tex"
)
```

## Save Tex File - Subset of Coefficients (Table 2)

```{r model_zc_interactions_white_college_income_kable}
model_zc_interactions_white_college_income_kable_COEF_SUBSET_fe_info <- 
  model_zc_interactions_white_college_income_kable_ready %>% 
  filter(display_name %in% c(
    "Republican",
    "Share GOP",
    "Republican * Share GOP",
    "Republican * Share White ('18 Zip)",
    "Republican * Share College and Above ('18 ZIP)",
    "Republican * Median HH Income, Thousands ('18 ZIP)",
    "N", "R-Squared", "Adj. R-Squared"
  )) %>% 
  create_kable_felm(
    data = .,
    dv_name_vec = dv_name_vec_four_primary,
    input_term_names = df_display_names_regressions, 
    packed_rows = FALSE, 
    escape_first_space = TRUE,
    caption = "Within-ZIP Code Models of Neighborhood Partisanship",
    fixed_effect_info = list(
      vars = names(dv_name_vec_four_primary),
      fe   = c("\\ Individual Controls", 
               "\\ Survey Wave Fixed Effect", 
               "\\ ZIP Code Fixed Effect"),
      cell_value = "\\ \\ \\ \\ \\ \\ \\checkmark"
    )
  ) %>% 
  footnote(general = "See Table S8 for individual controls and wave fixed effects",
           general_title =  "Full model:",
           footnote_as_chunk = T,
           title_format = "bold")
```

```{r SAVE_model_zc_interactions_white_college_income_kable_COEF_SUBSET_fe_info}
if (params$save_tex) {
  save_kable(
    x = model_zc_interactions_white_college_income_kable_COEF_SUBSET_fe_info,
    file = "tex-files/model_zc_interactions_white_college_income_kable_COEF_SUBSET_fe_info.tex"
  )
}
```

# White Interactions

## Run Model

```{r formula_zc_interactions_with_white_interaction}
formula_zc_interactions_with_white_interaction <- 
  dv_value ~ republican + 
             white:zip_white +
             white:prop_r +
             white:zip_educ_2018_college_and_above +
             white:zip_hh_inc_2018_median_range_middle +
             log(deaths_14_days_county_pc_1k + 0.1) +
             # Time Fixed Effects
             as.factor(wave) +
             # Individual Controls
             factor(demo_gender) +
             demo_household_income_range_middle +
             demo_household_income_missing +
             as_factor(demo_education_binned) +
             demo_age +
             # demo_race_ethnicity_binned_hispanic +
             white +
             Q7_health_dx_sum | 
             # FE + SE
             demo_zip_scrambled | 0 | demo_zip_scrambled
```

```{r model_zc_interactions_with_white_interaction}
model_zc_interactions_with_white_interaction <- data_for_model_list %>%
  # Run models
  split(.$dv_name) %>% 
  map( ~ felm(data = drop_na(., party),
              weights = drop_na(., party)$weights,
              formula = formula_zc_interactions_with_white_interaction))
```

## Save Tex File - All Coefficients (Table S9)

```{r model_zc_interactions_with_white_interaction_kable_fe_info}
model_zc_interactions_with_white_interaction_kable_ready <- 
  tidy_felm_prep_kable(
    input_model_list = model_zc_interactions_with_white_interaction,
    input_term_names =  df_display_names_regressions,
    dv_name_vec = dv_name_vec_four_primary
  )

model_zc_interactions_with_white_interaction_kable_fe_info <- create_kable_felm(
  data = model_zc_interactions_with_white_interaction_kable_ready,
  dv_name_vec = dv_name_vec_four_primary,
  input_term_names = df_display_names_regressions,
  packed_rows = FALSE, 
  escape_first_space = TRUE,
  caption = "Interactions Between White and ZIP Code-level Measures",
  fixed_effect_info = list(
    vars = names(dv_name_vec_four_primary),
    fe   = c("\\ ZIP Code Fixed Effect"),
    cell_value = "\\ \\ \\ \\ \\ \\ \\checkmark"
  )
)
```

```{r}
save_kable(
  x = model_zc_interactions_with_white_interaction_kable_fe_info,
  file = "tex-files/model_zc_interactions_with_white_interaction_kable_fe_info.tex"
)
```

\pagebreak

# College and Above Interactions

## Run Model

```{r formula_zc_interactions_with_college_interaction}
formula_zc_interactions_with_college_interaction <- 
  dv_value ~ republican + 
             college_and_above:zip_educ_2018_college_and_above +
             college_and_above:zip_white +
             college_and_above:prop_r +
             college_and_above:zip_hh_inc_2018_median_range_middle +
             log(deaths_14_days_county_pc_1k + 0.1) +
             # Time Fixed Effects
             as.factor(wave) +
             # Individual Controls
             factor(demo_gender) +
             demo_household_income_range_middle +
             demo_household_income_missing +
             # as_factor(demo_education_binned) +
             college_and_above +
             demo_age +
             demo_race_ethnicity_binned_hispanic +
             Q7_health_dx_sum | 
             # FE + SE
             demo_zip_scrambled | 0 | demo_zip_scrambled
```

```{r model_zc_interactions_with_college_interaction}
model_zc_interactions_with_college_interaction <- data_for_model_list %>%
  # Run models
  split(.$dv_name) %>% 
  map( ~ felm(data = drop_na(., party),
              weights = drop_na(., party)$weights,
              formula = formula_zc_interactions_with_college_interaction))
```

## Save Tex File - All Coefficients (Table S10)

```{r model5_basic_model_with_college_interaction_kable_fe_info}
model_zc_interactions_with_college_interaction_kable_ready <- 
  tidy_felm_prep_kable(
    input_model_list = model_zc_interactions_with_college_interaction,
    input_term_names =  df_display_names_regressions,
    dv_name_vec = dv_name_vec_four_primary
  )

model_zc_interactions_with_college_interaction_kable_fe_info <- create_kable_felm(
  data = model_zc_interactions_with_college_interaction_kable_ready,
  dv_name_vec = dv_name_vec_four_primary,
  input_term_names = df_display_names_regressions,
  packed_rows = FALSE, 
  escape_first_space = TRUE,
  caption = " Interactions Between College and Above and ZIP Code-level Measures",
  fixed_effect_info = list(
    vars = names(dv_name_vec_four_primary),
    fe   = c("\\ ZIP Code Fixed Effect"),
    cell_value = "\\ \\ \\ \\ \\ \\ \\checkmark"
  )
)
```

```{r}
save_kable(
  x = model_zc_interactions_with_college_interaction_kable_fe_info,
  file = "tex-files/model_zc_interactions_with_college_interaction_kable_fe_info.tex"
)
```


\pagebreak

# High-Income Interactions

## Run Model

```{r formula_zc_interactions_with_income_interaction}
formula_zc_interactions_with_income_interaction <- 
  dv_value ~ republican + 
             demo_household_income_top_quartile_binary:zip_hh_inc_2018_prop_top_quartile + 
             demo_household_income_top_quartile_binary:zip_educ_2018_college_and_above +
             demo_household_income_top_quartile_binary:zip_white +
             demo_household_income_top_quartile_binary:prop_r +
             log(deaths_14_days_county_pc_1k + 0.1) +
             # Time Fixed Effects
             as.factor(wave) +
             # Individual Controls
             factor(demo_gender) +
             # demo_household_income_range_middle +
             demo_household_income_top_quartile_binary +
             demo_household_income_missing +
             as_factor(demo_education_binned) +
             demo_age +
             demo_race_ethnicity_binned_hispanic +
             Q7_health_dx_sum | 
             # FE + SE
             demo_zip_scrambled | 0 | demo_zip_scrambled
```

```{r model_zc_interactions_with_income_interaction}
model_zc_interactions_with_income_interaction <- data_for_model_list %>%
  # Run models
  split(.$dv_name) %>% 
  map( ~ felm(data = drop_na(., party),
              weights = drop_na(., party)$weights,
              formula = formula_zc_interactions_with_income_interaction))
```

## Save Tex File - All Coefficients (Table S11)

```{r model_zc_interactions_with_income_interaction_kable_fe_info}
model_zc_interactions_with_income_interaction_kable_ready <- 
  tidy_felm_prep_kable(
    input_model_list = model_zc_interactions_with_income_interaction,
    input_term_names =  df_display_names_regressions,
    dv_name_vec = dv_name_vec_four_primary
  )

model_zc_interactions_with_income_interaction_kable_fe_info <- create_kable_felm(
  data = model_zc_interactions_with_income_interaction_kable_ready,
  dv_name_vec = dv_name_vec_four_primary,
  input_term_names = df_display_names_regressions,
  packed_rows = FALSE, 
  escape_first_space = TRUE,
  caption = "Interactions Between High Income and ZIP Code-level Measures",
  fixed_effect_info = list(
    vars = names(dv_name_vec_four_primary),
    fe   = c("\\ ZIP Code Fixed Effect"),
    cell_value = "\\ \\ \\ \\ \\ \\ \\checkmark"
  )
)
```

```{r}
save_kable(
  x = model_zc_interactions_with_income_interaction_kable_fe_info,
  file = "tex-files/model_zc_interactions_with_income_interaction_kable_fe_info.tex"
)
```

# Within-County Models

## Pooled Model (GOP * All Share Variables)

### Run Model

```{r formula_county_interactions_white_college_income}
formula_county_interactions_white_college_income <- 
  dv_value ~ republican + 
             prop_r + 
             republican:prop_r +
             republican:zip_white +
             republican:zip_educ_2018_college_and_above +
             republican:zip_hh_inc_2018_median_range_middle +
             log(deaths_14_days_county_pc_1k + 0.1) + 
             # Time fixed effects
             as.factor(wave) + 
             # Individual Controls
             factor(demo_gender) + 
             demo_household_income_range_middle + 
             demo_household_income_missing +      
             as_factor(demo_education_binned) + 
             demo_age + 
             demo_race_ethnicity_binned_hispanic + 
             Q7_health_dx_sum + 
             # Zipcode controls
             zip_pop_density_2018_over_2010 + 
             zip_educ_2018_college_and_above + 
             zip_white + 
             zip_hh_inc_2018_median_range_middle |   
             # FE + SE
             county_modified | 0 | demo_zip_scrambled
```

```{r model_county_interactions_white_college_income}
model_county_interactions_white_college_income <- data_for_model_list %>% 
  filter(!is.na(party)) %>% 
  # Run Models
  split(.$dv_name) %>% 
  map( ~ felm(data = .,
              weights = .$weights,
              formula = formula_county_interactions_white_college_income))
```

### Save Tex File - All Coefficients (Table S7)

```{r model_county_interactions_white_college_income_kable}
model_county_interactions_white_college_income_kable_ready <- 
  tidy_felm_prep_kable(
    input_model_list = model_county_interactions_white_college_income,
    input_term_names = df_display_names_regressions,
    dv_name_vec = dv_name_vec_four_primary
  )
```

```{r model_county_interactions_white_college_income_kable_fe_info}
model_county_interactions_white_college_income_kable_fe_info <-
  model_county_interactions_white_college_income_kable_ready %>% 
  create_kable_felm(
    data = .,
    dv_name_vec = dv_name_vec_four_primary,
    input_term_names = df_display_names_regressions,
    packed_rows = FALSE,
    escape_first_space = TRUE,
    caption = "Within-County Models of Neighborhood Partisanship (All Coefficients)",
    fixed_effect_info = list(
      vars = names(dv_name_vec_four_primary),
      fe   = c("\\ County Fixed Effect"),
      cell_value = "\\ \\ \\ \\ \\ \\ \\checkmark"
    )
  )
```

```{r}
if (params$save_tex) {
  save_kable(
    x = model_county_interactions_white_college_income_kable_fe_info,
    file = "tex-files/model_county_interactions_white_college_income_kable_fe_info.tex"
  )
}
```


\pagebreak

### Save Tex File - Subset of Coefficients (Table 1)

```{r model_county_interactions_white_college_income_kable_COEF_SUBSET}
model_county_interactions_white_college_income_kable_COEF_SUBSET <-
  model_county_interactions_white_college_income_kable_ready %>% 
  filter(display_name %in% c(
    "Republican",
    "Share GOP",
    "Republican * Share GOP",
    "Republican * Share White ('18 Zip)",
    "Republican * Share College and Above ('18 ZIP)",
    "Republican * Median HH Income, Thousands ('18 ZIP)",
    "N", "R-Squared", "Adj. R-Squared"
  )) %>% 
  create_kable_felm(
    data = .,
    dv_name_vec = dv_name_vec_four_primary,
    input_term_names = df_display_names_regressions,
    packed_rows = FALSE,
    escape_first_space = TRUE,
    caption = "Within-County Models of Neighborhood Partisanship",
    fixed_effect_info = list(
      vars = names(dv_name_vec_four_primary),
      fe   = c("\\ Individual Controls", 
               "\\ ZIP Code Controls",
               "\\ Survey Wave Fixed Effect", 
               "\\ County Fixed Effect"),
      cell_value = "\\ \\ \\ \\ \\ \\ \\checkmark"
    )
  ) %>% 
  footnote(general = "See Table S7 for individual controls and wave fixed effects",
           general_title =  "Full model:",
           footnote_as_chunk = T,
           title_format = "bold")
```

```{r}
if (params$save_tex) {
  save_kable(
    x = model_county_interactions_white_college_income_kable_COEF_SUBSET,
    file = "tex-files/model_county_interactions_white_college_income_kable_COEF_SUBSET.tex"
  )
}
```

\pagebreak

### Statistical Testing

```{r df_ci_prop_r_dem_gop}
# PULL OUT VARIANCE-COVARIANCE MATRIX FOR SHARE GOP AND THE INTERACTION TERM
model_county_interactions_white_college_income_vcov_party <- 
  model_county_interactions_white_college_income %>% 
  # Which dependent variable? (Mask-wearing)
  .[["Q31_mask_wearing_always_binary"]] %>% 
  # Find variance-covariance matrix
  vcov() %>% 
  # Pull out just the cells we need
  .[c("prop_r", "republican:prop_r"), 
    c("prop_r", "republican:prop_r")]

# CLEAN UP COEFFICIENTS AND "CROSS" WITH VALUES OF THE SHARE VARIABLE
df_ci_prop_r_republican_interaction_prep <- 
  model_county_interactions_white_college_income %>% 
  # Pull out and tidy mask model
  .[["Q31_mask_wearing_always_binary"]] %>% 
  tidy() %>% 
  # Filter
  filter(term %in% c("prop_r", "republican:prop_r")) %>% 
  # Pull out estimate
  select(term, estimate) %>% 
  spread(term, estimate) %>% 
  mutate(share_var_coef = prop_r,
         interaction_var_coef = `republican:prop_r`,
         model = "Share GOP") %>% 
  # Expand the dataset so each row is a value of the share variable
  crossing(
    tibble(share_var = 1)
  ) 

# CALCULATE EFFECT OF `prop.r` FOR REPUBLICANS
df_ci_prop_r_republican_interaction <- df_ci_prop_r_republican_interaction_prep %>% 
  mutate(share_plus_interaction_times_share = 
           share_var_coef * share_var + interaction_var_coef * share_var) %>%
  # Calculate standard error. Note: GOP = 1 in this calculation
  # x^2 * V(B1hat) + x^2 * GOP * V(B2hat) + 2 * x * GOP * Cov(B1hat, B2hat) 
  # To find SE: take square root
  mutate(Var_B1hat = model_county_interactions_white_college_income_vcov_party["prop_r", "prop_r"],
         Var_B2hat = model_county_interactions_white_college_income_vcov_party["republican:prop_r", "republican:prop_r"],
         Cov_B1hat_B2hat = model_county_interactions_white_college_income_vcov_party["prop_r", "republican:prop_r"],
         # Var_manual_by_components = component1 + component2 + component3,
         Var_manual = share_var^2 * Var_B1hat + 
                      share_var^2 * Var_B2hat + 
                      2 * share_var^2 * Cov_B1hat_B2hat, # Square share_var here!
         SE_manual = sqrt(Var_manual)) %>% 
  # Create confidence intervals
  mutate(party = "Republican")

# CALCULATE EFFECT OF `prop.r` FOR DEMOCRATS
df_ci_prop_r_alone <- df_ci_prop_r_republican_interaction_prep %>% 
  mutate(share_plus_interaction_times_share = share_var_coef * share_var) %>% 
  # Calculate standard error. Except now GOP = 0... so non-first terms fall out.
  # x^2 * V(B1hat) + x^2 * GOP * V(B2hat) + 2 * x * GOP * Cov(B1hat, B2hat) 
  # To find SE: take square root
  mutate(Var_B1hat = model_county_interactions_white_college_income_vcov_party["prop_r", "prop_r"],
         Var_manual = share_var^2 * Var_B1hat, 
         SE_manual = sqrt(Var_manual)) %>% 
  mutate(party = "Democrat")

# COMBINE DATA FOR PLOT AND CALCULATE CONFIDENCE INTERVALS
df_ci_prop_r_dem_gop <- bind_rows(
  df_ci_prop_r_alone,
  df_ci_prop_r_republican_interaction
) %>% 
  mutate(ci.low = share_plus_interaction_times_share + qnorm(0.025) * SE_manual,
         ci.high = share_plus_interaction_times_share + qnorm(0.975) * SE_manual)
```

```{r g_ci_prop_r_dem_gop, fig.width=4, fig.height=4}
g_ci_prop_r_dem_gop_single_point <- df_ci_prop_r_dem_gop %>% 
  mutate(party = fct_relevel(party, "Republican")) %>% 
  ggplot(aes(x = party, y = share_plus_interaction_times_share)) +
  geom_hline(yintercept = 0,
             color = "darkgray") +
  geom_point(aes(color = party,
                fill = party)) +
  geom_errorbar(aes(ymin = ci.low,
                  ymax = ci.high,
                  color = party,
                  fill = party),
                width = 0.2,
              ) +
  # Theme and Structure
  theme_minimal() +
  theme(legend.position = "none") +
  # Scales and Labels
  labs(x = NULL, y = "Effect of Share GOP on Wearing a Mask") +
  expand_limits(y = c(-0.4, 0.4)) +
  scale_color_manual(values = colors, name = NULL) +
  scale_fill_manual(values = colors, name = NULL) 

g_ci_prop_r_dem_gop_single_point
```

```{r}
ggsave(
  plot = g_ci_prop_r_dem_gop_single_point,
  filename = "Figures/Figure 3.pdf",
  width = 11,
  height = 11, 
  units = "cm"
)
```

# Figure 4: Within-Zipcode Rate Plot

```{r}
df_zip_party_dv <- dat_long_dvs %>% 
  filter(party_with_independents != "Independent",
         dv_name == "mask_wearing_always") %>% 
  group_by(dv_name, demo_zip_scrambled, prop_r, party_with_independents) %>% 
  summarise(w_prop = weighted.mean(x = dv_value, 
                                   w = weights,
                                   na.rm = TRUE),
            .groups = "drop") %>% 
  spread(party_with_independents, w_prop) %>% 
  mutate(dem_minus_rep = Democrat - Republican,
         rep_minus_dem = Republican - Democrat)
```

## Republican - Democrat with Modeled Regression Line

```{r}
df_reg_line_figure4_rep_minus_dem <- tibble(
    prop_r = seq(0, 1, 1),
    prop_r2 = seq(0, 1, 1)
  ) %>% 
  mutate(rep_minus_dem = prop_r * model_zc_interactions_white_college_income %>% 
                .$Q31_mask_wearing_always_binary %>% 
                coef() %>% 
                .["republican:prop_r"]) %>% 
  mutate(rep_minus_dem = rep_minus_dem + model_zc_interactions_white_college_income %>% 
                .$Q31_mask_wearing_always_binary %>% 
                coef() %>% 
                .["republican"]) %>% 
  gather(key, val, - prop_r2) %>% 
  unite(key, key, prop_r2) %>% 
  spread(key, val)
```

```{r}
g_zip_party_dv_rep_minus_dem <- df_zip_party_dv %>% 
  drop_na(dem_minus_rep) %>% 
  # Plot
  ggplot(aes(x = prop_r, y = rep_minus_dem)) +
  geom_hline(yintercept = 0, color = "darkgray") +
  stat_summary_bin(aes(x = prop_r, y = rep_minus_dem),
                   fun.data = mean_ci,
                   breaks = seq(0, 0.9, 0.1),
                   color = "darkgray") +
  geom_smooth(method = "loess", alpha = 0.2, size = 0.5, color = "darkgray") + 
  geom_segment(data = df_reg_line_figure4_rep_minus_dem,
               aes(x = prop_r_0,
                   xend = prop_r_1,
                   y = rep_minus_dem_0,
                   yend = rep_minus_dem_1),
               color = "black",
               size = 1,
               linetype = "dashed") +
  # Theme and Structure
  theme_minimal() +
  theme(legend.position = "bottom") +
  # Labels and Scales
  labs(y = "Percentage Point Difference\n(Republican - Democrat) ",
       x = '% GOP in ZIP Code') +
  scale_y_continuous(labels = percent) + #, breaks = seq(-1, 1, 0.2)) +
  scale_x_continuous(labels = percent, expand = expand_scale(mult = 0.1)) +
  expand_limits(y = c(-0.3, 0.3), x = 1) +
  scale_color_manual(values = colors, name = NULL) +
  scale_fill_manual(values = colors, name = NULL)

g_zip_party_dv_rep_minus_dem
```

```{r}
ggsave(
  filename = "Figures/Figure 4.pdf",
  plot = g_zip_party_dv_rep_minus_dem,
  width = 11,  
  height = 11, 
  units = "cm",
  dpi = 300
)
```


# Supplement Results: Ideology (Tables S12 - S13)

## Within-ZIP Code

### Wave 4 Only, No Ideology

```{r formula_ideo1_wave4_no_ideology}
formula_ideo1_wave4_no_ideology <- 
  dv_value ~ republican + 
             # prop_r + 
             republican:prop_r +
             republican:zip_white +
             republican:zip_educ_2018_college_and_above +
             republican:zip_hh_inc_2018_median_range_middle +
             log(deaths_14_days_county_pc_1k + 0.1) +
             # Time fixed effects not needed (1 wave only)
             # Individual Controls
             factor(demo_gender) + 
             demo_household_income_range_middle + 
             demo_household_income_missing +      
             as_factor(demo_education_binned) + 
             demo_age + 
             demo_race_ethnicity_binned_hispanic + 
             Q7_health_dx_sum | 
             # FE + SE
             demo_zip_scrambled | 0 | demo_zip_scrambled
```

Models not shown (see later table)

```{r model_ideo1_wave4_no_ideology}
model_ideo1_wave4_no_ideology <- data_for_model_list %>% 
  filter(wave == 4) %>% 
  filter(!is.na(party)) %>%
  # Run Models
  split(.$dv_name) %>% 
  map( ~ felm(data = .,
              weights = .$weights,
              formula = formula_ideo1_wave4_no_ideology))
```

### Wave 4 Only, Ideology Included

```{r formula_ideo2_wave4_ideology_included}
formula_ideo2_wave4_ideology_included <- 
  dv_value ~ republican + 
             # prop_r + 
             republican:prop_r +
             republican:zip_white +
             republican:zip_educ_2018_college_and_above +
             republican:zip_hh_inc_2018_median_range_middle +
             log(deaths_14_days_county_pc_1k + 0.1) +
             # Ideology
             demo_ideology_f + 
             # Time fixed effects not needed (1 wave only)
             # Individual Controls
             factor(demo_gender) + 
             # demo_household_income_no_unanswered +  
             demo_household_income_range_middle + 
             demo_household_income_missing +      
             as_factor(demo_education_binned) + 
             demo_age + 
             demo_race_ethnicity_binned_hispanic + 
             Q7_health_dx_sum |  
             # FE + SE
             demo_zip_scrambled | 0 | demo_zip_scrambled
```

Models not shown (see later table)

```{r model_ideo2_wave4_ideology_included}
model_ideo2_wave4_ideology_included <- data_for_model_list %>% 
  filter(wave == 4) %>% 
  filter(!is.na(party)) %>%
  # Run Models
  split(.$dv_name) %>% 
  map( ~ felm(data = .,
              weights = .$weights,
              formula = formula_ideo2_wave4_ideology_included))
```

### Save Tex-File - All Coefficients

```{r}
model_ideo1_ideo2_masks_covid_vaccines <- list(
  masks_no_ideology   = model_ideo1_wave4_no_ideology$Q31_mask_wearing_always_binary,
  masks_with_ideology = model_ideo2_wave4_ideology_included$Q31_mask_wearing_always_binary,
  covid_vaccine_no_ideology   = model_ideo1_wave4_no_ideology$Q147_would_get_covid_vaccine_yes_binary,
  covid_vaccine_with_ideology = model_ideo2_wave4_ideology_included$Q147_would_get_covid_vaccine_yes_binary
)
```


```{r}
dv_name_vec_four_primary_w_wo_ideology <- c(
  "masks_no_ideology"           = "Without Ideology",
  "masks_with_ideology"         = "With Ideology",
  "covid_vaccine_no_ideology"   = "Without Ideology",
  "covid_vaccine_with_ideology" = "With Ideology"
)

model_ideo1_ideo2_masks_covid_vaccines_kable_ready <- 
  tidy_felm_prep_kable(
    input_model_list = model_ideo1_ideo2_masks_covid_vaccines,
    input_term_names = df_display_names_regressions,
    dv_name_vec = dv_name_vec_four_primary_w_wo_ideology
  ) 


model_ideo1_ideo2_masks_covid_vaccines_kable <-
  model_ideo1_ideo2_masks_covid_vaccines_kable_ready %>% 
  create_kable_felm(
    data = .,
    dv_name_vec = dv_name_vec_four_primary_w_wo_ideology,
    input_term_names = df_display_names_regressions,
    packed_rows = FALSE,
    escape_first_space = TRUE,
    caption = "Comparing Within-ZIP Code Models With and Without Ideology (December 4-16 Wave Only)",
    fixed_effect_info = list(
      vars = names(dv_name_vec_four_primary_w_wo_ideology),
      fe   = c("\\ ZIP Code Fixed Effect"),
      cell_value = "\\ \\ \\ \\ \\ \\ \\checkmark"
    )
  ) %>% 
  add_header_above(c(" " = 1, 
                     "Always Wear Mask" = 2,
                     "COVID Vaccine" = 2))
```

```{r}
if (params$save_tex) {
  save_kable(
    x = model_ideo1_ideo2_masks_covid_vaccines_kable,
    file = "tex-files/model_ideo1_ideo2_masks_covid_vaccines_kable.tex"
  )
}
```

## Within-County

### Wave 4 Only, No Ideology

```{r formula_ideo1_wave4_no_ideology_COUNTY}
formula_ideo1_wave4_no_ideology_COUNTY <- 
  dv_value ~ republican + 
             prop_r + 
             republican:prop_r +
             republican:zip_white +
             republican:zip_educ_2018_college_and_above +
             republican:zip_hh_inc_2018_median_range_middle +
             log(deaths_14_days_county_pc_1k + 0.1) + 
             # Time fixed effects not needed (1 wave only)
             # Individual Controls
             factor(demo_gender) + 
             demo_household_income_range_middle + 
             demo_household_income_missing +      
             as_factor(demo_education_binned) + 
             demo_age + 
             demo_race_ethnicity_binned_hispanic + 
             Q7_health_dx_sum + 
             # Zipcode controls
             zip_pop_density_2018_over_2010 + 
             zip_educ_2018_college_and_above + 
             zip_white + 
             zip_hh_inc_2018_median_range_middle |   
             # FE + SE
             county_modified | 0 | demo_zip_scrambled
```

```{r model_set_14_pooled_gop_all_share_county_list_county_fe_zipcode_se}
model_ideo1_wave4_no_ideology_COUNTY <- data_for_model_list %>% 
  filter(wave == 4) %>% 
  filter(!is.na(party)) %>% 
  # Run Models
  split(.$dv_name) %>% 
  map( ~ felm(data = .,
              weights = .$weights,
              formula = formula_ideo1_wave4_no_ideology_COUNTY))
```

### Wave 4 Only, Ideology Included

```{r formula_ideo2_wave4_ideology_included_COUNTY}
formula_ideo2_wave4_ideology_included_COUNTY <- 
  dv_value ~ republican + 
             prop_r + 
             republican:prop_r +
             republican:zip_white +
             republican:zip_educ_2018_college_and_above +
             republican:zip_hh_inc_2018_median_range_middle +
             log(deaths_14_days_county_pc_1k + 0.1) + 
             # Ideology
             demo_ideology_f + 
             # Time fixed effects not needed (1 wave only)
             factor(demo_gender) + 
             demo_household_income_range_middle +
             demo_household_income_missing +     
             as_factor(demo_education_binned) + 
             demo_age + 
             demo_race_ethnicity_binned_hispanic + 
             Q7_health_dx_sum +
             # Zipcode controls
             zip_pop_density_2018_over_2010 + 
             zip_educ_2018_college_and_above + 
             zip_white + 
             zip_hh_inc_2018_median_range_middle |    
             # FE + SE
             county_modified | 0 | demo_zip_scrambled
```

```{r model_ideo2_wave4_ideology_included_COUNTY}
model_ideo2_wave4_ideology_included_COUNTY <- data_for_model_list %>% 
  filter(wave == 4) %>% 
  filter(!is.na(party)) %>% 
  # Run Models
  split(.$dv_name) %>% 
  map( ~ felm(data = .,
              weights = .$weights,
              formula = formula_ideo2_wave4_ideology_included_COUNTY))
```

### Save Tex Files - All Coefficients

```{r}
model_ideo1_ideo2_masks_covid_vaccines_COUNTY <- list(
  masks_no_ideology   = model_ideo1_wave4_no_ideology_COUNTY$Q31_mask_wearing_always_binary,
  masks_with_ideology = model_ideo2_wave4_ideology_included_COUNTY$Q31_mask_wearing_always_binary,
  covid_vaccine_no_ideology   = model_ideo1_wave4_no_ideology_COUNTY$Q147_would_get_covid_vaccine_yes_binary,
  covid_vaccine_with_ideology = model_ideo2_wave4_ideology_included_COUNTY$Q147_would_get_covid_vaccine_yes_binary
)
```


```{r}
model_ideo1_ideo2_masks_covid_vaccines_kable_ready_COUNTY <- 
  tidy_felm_prep_kable(
    input_model_list = model_ideo1_ideo2_masks_covid_vaccines_COUNTY,
    input_term_names = df_display_names_regressions,
    dv_name_vec = dv_name_vec_four_primary_w_wo_ideology
  ) 


model_ideo1_ideo2_masks_covid_vaccines_kable_COUNTY <-
  model_ideo1_ideo2_masks_covid_vaccines_kable_ready_COUNTY  %>% 
  create_kable_felm(
    data = .,
    dv_name_vec = dv_name_vec_four_primary_w_wo_ideology,
    input_term_names = df_display_names_regressions,
    packed_rows = FALSE,
    escape_first_space = TRUE,
    caption = "Comparing Within-County Models With and Without Ideology (December 4-16 Wave Only)",
    fixed_effect_info = list(
      vars = names(dv_name_vec_four_primary_w_wo_ideology),
      fe   = c("\\ County Fixed Effect"),
      cell_value = "\\ \\ \\ \\ \\ \\ \\checkmark"
    )
  ) %>% 
  add_header_above(c(" " = 1, 
                     "Always Wear Mask" = 2,
                     "COVID Vaccine" = 2))
```

```{r}
if (params$save_tex) {
  save_kable(
    x = model_ideo1_ideo2_masks_covid_vaccines_kable_COUNTY,
    file = "tex-files/model_ideo1_ideo2_masks_covid_vaccines_kable_COUNTY.tex"
  )
}
```

